O’Hara on GitHub


knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
                      fig.height = 4, fig.width = 7)

library(terra)
library(sf)
library(oharac)
library(data.table)
library(tidyverse)
library(here)
library(cowplot)
source(here('common_fxns.R'))

1 Summary

Sensitivity analysis to determine how drastically scores change with different values of functional traits.

In script 4a, we calculated sensitivity around functional diversity metrics and per-stressor impacts. In 4b, we aggregate results across all stressors for each randomized iteration, to determine the distribution of mean values for cumulative impacts. Here, we map them out and make cool figures.

2 Methods

For each randomized trait, plot the variation in functional diversity metrics for the sampled points, with a nearest-neighbor interpolation to help spot regional variation. We will examine several stats: absolute and relative difference (diff / mean), 95% confidence interval, coefficient of variance (std dev divided by mean)

ocean_r <- rast(here('_spatial/ocean_area_mol.tif'))
cell_xy <- ocean_r %>% 
  setValues(1:ncell(.)) %>%
  setNames('cell_id') %>%
  as.data.frame(xy = TRUE)
fd_fs <- list.files(here_anx('sens_analysis/funct_div_summary'), 
                    pattern = '.csv', full.names = TRUE)
fd_ts <- basename(fd_fs) %>% str_remove_all('.+_shuffle_|.csv')


fd_df <- lapply(fd_fs, fread) %>%
  setNames(fd_ts) %>%
  rbindlist(idcol = 't') %>%
  mutate(n_fe_diff = n_fe_mean - n_fe,
         fv_diff = fv_mean - fv,
         n_fe_ci = (n_fe_95hi - n_fe_95lo),
         fv_ci   = (fv_95hi - fv_95lo),
         n_fe_cv = n_fe_sdev / n_fe_mean,
         fv_cv   = fv_sdev / fv_mean) %>%
  rename(n_fe_value = n_fe, fv_value = fv) %>%
  pivot_longer(cols = -c(cell_id, t, nspp),
               names_to = 'var', values_to = 'val') %>%
  mutate(param = str_extract(var, 'n_fe|fv'),
         stat  = str_extract(var, 'mean|sdev|ci|cv|diff|95..|value')) %>%
  select(-var) %>%
  left_join(cell_xy, by = 'cell_id')

2.1 Functional diversity metrics

2.1.1 Maps of functional vulnerability

How does functional vulnerability vary as traits are randomized?

land_sf <- read_sf(here('_spatial/ne_10m_land/ne_10m_land_no_casp.shp')) %>%
  st_transform(st_crs(ocean_r))

globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                      c(180, 90), c(180, -90), c(-180, -90)) 
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = crs(ocean_r))
map_param <- function(df, s, p, r = ocean_r, sig = FALSE, thresh = NULL) {

  val_mtx <- df %>%
    filter(stat == s & param == p) %>%
    select(x, y, z = val) %>%
    as.matrix()
  
  val_rast <- interpNear(r, val_mtx, radius = 500000) %>%
    mask(r)
  
  if(sig) {
    sig_mtx <- df %>%
      filter(param == p & stat %in% c('95hi', '95lo', 'value')) %>%
      pivot_wider(names_from = stat, values_from = val) %>%
      mutate(sig = ifelse(between(value, `95lo`, `95hi`), 0, 1)) %>%
      select(x, y, z = sig) %>%
      as.matrix()
    sig_rast <- interpNear(r, sig_mtx, radius = 50000) %>%
      mask(r)
    val_rast[sig_rast == 0] <- NA
  }
  if(!is.null(thresh)) val_rast[abs(val_rast) < thresh] <- NA

  return(val_rast)
}

finish_map <- function(p, div = TRUE) {
  ### div is whether the palette should be a diverging palette (TRUE) or
  ### sequential (FALSE)

  p <- p +
    ### continents:
    geom_sf(data = land_sf,
            fill = 'grey96', color = 'grey40',
            size = .10) +
    geom_sf(data = globe_border,
            fill = NA, color = 'grey70', 
            size = .1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    facet_wrap(~ trait) +
    theme_void()
  
  if(div) {
    ### diverging palette
    p <- p + 
      scale_fill_gradient2(low = '#d01c8b', mid = '#ffffdf', high = '#4dac26')
  } else {
    ### sequential palette
    p <- p + scale_fill_viridis_c()
  }
  
  return(p)
}
fv_fig_f <- '_figs/fv_diff_by_trait_map.png'

if(!file.exists(fv_fig_f)) {
  fv_map_list <- vector('list', length = length(fd_ts)) %>%
    setNames(fd_ts)
  for(trait in fd_ts) {
    # trait <- fd_ts[1]
    fv_map_list[[trait]] <- map_param(df = fd_df %>% filter(t == trait),
                                      p = 'fv', s = 'diff')
  }
  
  fv_map_stack <- rast(fv_map_list)
  mean_sd_vals <- as.data.frame(fv_map_stack) %>%
    summarize(across(everything(), .fns = list(mean = mean, sd = sd), .names = '{.col}_{.fn}'))
  #     len_mean     len_sd   mob_mean     mob_sd   trp_mean     trp_sd  wcol_mean    wcol_sd
  # 1 0.03633012 0.04194202 0.03743112 0.04292291 0.05350717 0.04770103 0.02206068 0.02724206  
  
  fv_map_df <- as.data.frame(fv_map_stack, xy = TRUE) %>%
    pivot_longer(-c(x, y), names_to = 'trait', values_to = 'delta_fv') %>%
    mutate(trait = case_when(trait == 'len' ~ 'body length',
                             trait == 'mob' ~ 'adult mobility',
                             trait == 'trp' ~ 'trophic level',
                             trait == 'wcol' ~ 'water column position'))
    
  p <- ggplot() +
    geom_raster(data = fv_map_df, aes(x, y, fill = delta_fv)) +
    labs(fill = 'ΔFV')
  
  p <- finish_map(p)
  
  ggsave(fv_fig_f, width = 6.5, height = 3, dpi = 300)
  
}
knitr::include_graphics(fv_fig_f)

  # tbd
  ### param dist plot
  # tbd
  ### spp vs param plot
fv2_fig_f <- '_figs/fv_cv_by_trait_map.png'

if(!file.exists(fv2_fig_f)) {
  fv_map_list <- vector('list', length = length(fd_ts)) %>%
    setNames(fd_ts)
  for(trait in fd_ts) {
    fv_map_list[[trait]] <- map_param(df = fd_df %>% filter(t == trait), 
                                      p = 'fv', s = 'cv')
  }
  
  fv_map_stack <- rast(fv_map_list)
  
  fv_map_df <- as.data.frame(fv_map_stack, xy = TRUE) %>%
    pivot_longer(-c(x, y), names_to = 'trait', values_to = 'cv_fv') %>%
    mutate(trait = case_when(trait == 'len' ~ 'body length',
                             trait == 'mob' ~ 'adult mobility',
                             trait == 'trp' ~ 'trophic level',
                             trait == 'wcol' ~ 'water column position'))
    
  p <- ggplot() +
    geom_raster(data = fv_map_df, aes(x, y, fill = cv_fv)) +
    labs(fill = 'CV(FV)')
  
  p <- finish_map(p, div = FALSE)
  
  ggsave(fv2_fig_f, width = 6.5, height = 3, dpi = 300)

}
knitr::include_graphics(fv2_fig_f)

2.1.2 Plots of functional diversity parameters vs. number of species

How do functional diversity metrics (FV and number of FEs) vary as traits are randomized, relative to the number of species in a given pixel?

plot_f <- '_figs/fv_stat_by_trait_scatter.png'

if(!file.exists(plot_f)) {
  fv_plot_df <- fd_df %>%
    filter(stat %in% c('mean', 'cv', 'diff')) %>%
    filter(param == 'fv') %>%
    mutate(metric = paste(stat, param)) %>%
    mutate(metric = str_replace(metric, 'fv', '(FV)'),
           metric = str_replace(metric, 'mean ', 'mean'),
           metric = str_replace(metric, 'cv ', 'CV'),
           metric = str_replace(metric, 'diff ', 'Δ')) %>%
    mutate(t = case_when(t == 'len' ~ 'body length',
                         t == 'mob' ~ 'adult mobility',
                         t == 'trp' ~ 'trophic level',
                         t == 'wcol' ~ 'water column position')) %>%
    sample_frac(.1)
  
  p <- ggplot(fv_plot_df, aes(x = nspp, y = val)) +
    geom_point(alpha = .2) +
    theme_ohara() +
    facet_grid(metric ~ t, scales = 'free_y') +
    labs(x = 'Species richness') +
    theme(axis.title.y = element_blank())

    ggsave(plot_f, width = 6, height = 3, dpi = 300)
}

knitr::include_graphics(plot_f)

plot_f <- '_figs/n_fe_stat_by_trait_scatter.png'

if(!file.exists(plot_f)) {
  
  fe_plot_df <- fd_df %>%
    filter(stat %in% c('mean', 'cv', 'diff')) %>%
    filter(param == 'n_fe') %>%
    mutate(metric = paste(stat, param)) %>%
    mutate(metric = str_replace(metric, 'n_fe', '(# FE)'),
           metric = str_replace(metric, 'mean ', 'mean'),
           metric = str_replace(metric, 'cv ', 'CV'),
           metric = str_replace(metric, 'diff ', 'Δ')) %>%
    mutate(t = case_when(t == 'len' ~ 'body length',
                         t == 'mob' ~ 'adult mobility',
                         t == 'trp' ~ 'trophic level',
                         t == 'wcol' ~ 'water column position')) %>%
    sample_frac(.1)
  
  p <- ggplot(fe_plot_df, aes(x = nspp, y = val)) +
    geom_point(alpha = .2) +
    theme_ohara() +
    facet_grid(metric ~ t, scales = 'free_y') +
    labs(x = 'Species richness') +
    theme(axis.title.y = element_blank())
  
    ggsave(plot_f, width = 6, height = 3, dpi = 300)
}

knitr::include_graphics(plot_f)

2.2 Cumulative impact metrics

2.2.1 Maps of cumulative impact variation

How does cumulative impact vary as traits are randomized?

chi_fs <- list.files(here_anx('sens_analysis/impact_summary'), 
                    pattern = 'chi_.+.csv', full.names = TRUE)
chi_ts <- basename(chi_fs) %>% str_remove_all('.+_shuffle_|.csv')

chi_true_df <- rast(here('_output/cumulative_impact_maps/chi_funct_entity.tif')) %>%
  as.data.frame(xy = TRUE) %>%
  oharac::dt_join(cell_xy, by = c('x', 'y'), type = 'left') %>%
  filter(cell_id %in% fd_df$cell_id)

nspp_df <- fd_df %>%
  select(cell_id, nspp) %>%
  distinct()

chi_df <- lapply(chi_fs, fread) %>%
  setNames(chi_ts) %>%
  rbindlist(idcol = 't') %>%
  left_join(chi_true_df, by = 'cell_id') %>%
  left_join(nspp_df, by = 'cell_id') %>%
  mutate(chi_diff = chi_mean - chi_fe,
         chi_ci = (chi_95hi - chi_95lo),
         chi_cv = chi_sdev / chi_mean) %>%
  select(cell_id, t, nspp,
         chi_mean, chi_diff,
         chi_ci, chi_cv) %>%
  pivot_longer(cols = -c(cell_id, t, nspp),
               names_to = 'var', values_to = 'val') %>%
  mutate(param = 'chi',
         stat  = str_extract(var, 'mean|sdev|ci|cv|diff')) %>%
  select(-var) %>%
  left_join(cell_xy, by = 'cell_id')
chi_fig_f <- '_figs/chi_diff_by_trait_map.png'

if(!file.exists(chi_fig_f)) {
  chi_map_list <- vector('list', length = length(chi_ts)) %>%
    setNames(chi_ts)
  for(trait in fd_ts) {
    chi_map_list[[trait]] <- map_param(df = chi_df %>% filter(t == trait), 
                                       p = 'chi', s = 'diff')
  }
  
  chi_map_stack <- rast(chi_map_list)
  mean_sd_vals <- as.data.frame(chi_map_stack) %>%
    summarize(across(everything(), .fns = list(mean = mean, sd = sd), .names = '{.col}_{.fn}'))
  #      len_mean     len_sd     mob_mean     mob_sd    trp_mean     trp_sd    wcol_mean    wcol_sd
  # 1 -0.00801914 0.01990471 0.0008610911 0.01656777 0.003189031 0.01596556 -0.006297416 0.02280766  
  chi_map_df <- as.data.frame(chi_map_stack, xy = TRUE) %>%
    pivot_longer(-c(x, y), names_to = 'trait', values_to = 'delta_chi') %>%
    mutate(trait = case_when(trait == 'len' ~ 'body length',
                             trait == 'mob' ~ 'adult mobility',
                             trait == 'trp' ~ 'trophic level',
                             trait == 'wcol' ~ 'water column position'))
    
  p <- ggplot() +
    geom_raster(data = chi_map_df, aes(x, y, fill = delta_chi)) +
    labs(fill = 'ΔCHI')
  
  p <- finish_map(p)
  
  ggsave(chi_fig_f, width = 6.5, height = 3, dpi = 300)
  
}
knitr::include_graphics(chi_fig_f)

2.2.2 Plots of CHI stats vs. number of species

How do CHI stats vary as traits are randomized, relative to the number of species in a given pixel?

plot_f <- '_figs/chi_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
  chi_plot_df <- chi_df %>%
    filter(stat %in% c('mean', 'cv', 'diff')) %>%
    mutate(metric = paste(stat, param)) %>%
    mutate(metric = str_replace(metric, 'chi', '(CHI)'),
           metric = str_replace(metric, 'mean ', 'mean'),
           metric = str_replace(metric, 'cv ', 'CV'),
           metric = str_replace(metric, 'diff ', 'Δ')) %>%
    mutate(t = case_when(t == 'len' ~ 'body length',
                         t == 'mob' ~ 'adult mobility',
                         t == 'trp' ~ 'trophic level',
                         t == 'wcol' ~ 'water column position')) %>%
    sample_frac(.1)
  
  p <- ggplot(chi_plot_df, aes(x = nspp, y = val)) +
    geom_point(alpha = .2) +
    theme_ohara() +
    facet_grid(metric ~ t, scales = 'free_y') +
    labs(x = 'Species richness') +
    theme(axis.title.y = element_blank())
  
  ggsave(plot_f, width = 6, height = 3, dpi = 300)
}

knitr::include_graphics(plot_f)

3 Fig for SI

Combine the FV variation and CHI variation into a large panel.

p <- ggdraw() +
  draw_image(fv_fig_f,  x = 0.0, y = 0.51, width = 1.0, height = 0.49) +
  draw_image(chi_fig_f, x = 0.0, y = 0.0,  width = 1.0, height = 0.49) +
  draw_label('A.', x = .07, y = .995, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('B.', x = 0.5, y = .995, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('C.', x = .07, y = .748, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('D.', x = 0.5, y = .748, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('E.', x = .07, y = .485, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('F.', x = 0.5, y = .485, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('G.', x = .07, y = .235, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
  draw_label('H.', x = 0.5, y = .235, hjust = 0, vjust = 1,
             size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold')

si_fig_file <- here('5_ms_figs/figS6_sens_analysis.png')
ggsave(si_fig_file, height = 11.5, width = 12, units = 'cm', dpi = 300)

knitr::include_graphics(si_fig_file)